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INTRODUCTION 

A  series  of  three  recent  papers  (Snyder  and  Thomas,  1983;  Thomas 
and  Snyder,  1984a;  1984b)  has  introduced  methodology  for  performing 
probabilistic  analysis  and  synthesis  of  experimental  data  without 
assuming  that  a  particular  distribution  function  is  appropriate  to  a 
particular  data  set.  The  new  method  is  mathematical ly  form-free. 
Statistical  samples  are  smoothed  by  the  piece-wise  sliding  polynomials, 
yielding  form-free  probability  distributions.  These  probability  curves 
can  be  used  for  essentially  all  of  the  applications  purposes  of  conven¬ 
tional  distributions,  such  as  normal,  log-normal,  or  gamma  functions. 
The  methods  are  computer  oriented.  Therefore,  this  research  report 
has  been  prepared  for  the  convenience  of  potential  users.  General 
background  material  on  sliding  polynomials  may  be  found  in  earlier 
papers  by  Snyder,  1976  and  Snyder,  1980. 

Discussion  of  the  smoothing  methodology  will  be  brief,  intended 
only  to  provide  understanding  sufficient  for  use  and  adaptation  of 
programs.  The  user  is  referred  to  the  references  for  fuller  under¬ 
standing  . 

Smoothing  of  statistical  samples  with  sliding  polynomials  requires 
a  mathematical  transform  of  the  scale  of  the  original  data,  called 
variate  h,  to  a  scaled  abstract  variate  called  v.  This  transformation 
was  developed  so  that  the  smoothing  can  incorporate  boundary  controls 
at  infinite  values  of  h.  This  transformation  is  gradually  being  stan¬ 
dardized  through  additional  research.  The  programs  covered  in  this 
report  are  based  on  the  standardized  transform. 
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METHOD  OF  ANALYSIS 

The  method  of  smoothing  statistical  samples  with  sliding  polynomials 
and  with  incorporation  of  boundary  controls  is  shown  geometrically  in 
Figure  1.  Most  real  hydrologic  data  can  vary  between  zero  and  +  °°.  An 
exception  is  temperature,  where  real  values  may  be  negative.  With  the 
usual  zero  discontinuity,  the  0  to  +  °°  scale  for  the  hydrologic  variable, 
h,  is  shown  at  the  bottom  right  of  Figure  1.  h  is  first  converted  to 
a  standardized  variate,  h'  as  in  equation  1. 

h'  =  in 

ks 

in  is  the  mean  of  h,  s  is  the  standard  deviation  of  h,  and  k  is  an  em¬ 
pirical  parameter.  A  value  of  2  for  k  works  for  many  data  sets. 

h'  is  converted  to  v  by  equation  2. 

v  =  4.0  -  2.5  exp  (0.91 629h ' )  h  <  h 

h'  £  0 

v  >1.5  (2) 

v  =  1 .5  exp  (-1  . 5271  5h ' )  h  >  h 

h '  >_  0 
v  <1.5 

Equation  2  expresses  a  compound  curve  composed  of  two  descending  ex¬ 
ponentials  common  and  tangent  at  h'  =  0,  v  =  1.5.  For  calibration,  the 
curve  is  required  to  pass  through  the  point  h'  =  -1.0,  v  =  3.0. 

In  the  upper  section  of  Figure  1,  it  can  be  seen  that  three  arcs, 
or  spans,  of  the  sliding  polynomials  are  made  to  cover  the  possible 
range  of  h  from  zero  to  +  °°.  When  h  is  zero,  v  is  given  by  equation  3. 

v(h  =  0)  =  4  -  2.5  exp  (-0.91629  h/ks)  (3) 
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v  (h  =  0)  will  vary  from  sample  to  sample  as  h  and  s  vary.  The  interval 
from  v  =  0  to  v  (h  =  0)  is  divided  into  three  equal  segments  to  place 
nodes  at  uniform  interval  in  v-scale.  Nodes  0,  1,  2,  and  3  are  so 
located.  Using  the  same  interval,  boundary  nodes  are  placed  at  -1,  4, 
and  5.  The  location  of  these  boundary  nodes  has  no  meaning  in  v-scale, 
or  h-scale.  They  are  necessary  to  shape  the  sliding  polynomial  arcs. 

Class  limits  for  the  sample  are  defined  by  dividing  each  of  the 
three  polynomial  spans,  0-1,  1-2,  and  2-3,  into  tenths,  yielding  30 
classes  of  uniform  width  in  v-scale.  These  class  limits  then  are  con¬ 
verted  from  v-values  back  to  h-values  with  equations  1  and  2. 

The  statistical  sample  is  tallied  into  the  30  classes.  The  classes 
are  subtotalled,  starting  from  large  h  back  toward,  but  not  including, 
zero.  The  subtotals,  divided  by  the  sample  size,  and  multiplied  by 
100,  gives  a  sample  probability  of  any  individual  h  being  greater  than 
the  smaller  class  limit  in  h.  The  sample  is  then  ready  for  smoothing 
by  least  squares. 

Sliding  polynomial  smoothing  provides  values  of  the  nodes  at  1 ,  2, 
3,  and  4  on  the  node  placement  scale.  The  nodal  value  at  zero  place¬ 
ment  must  be  zero,  since  the  smoothing  curve  must  rise  from  zero  when 
h  is  +  An  additional  boundary  condition,  that  the  slope  of  the 
smoothing  curve  must  be  zero  when  h  is  +  is  imposed  on  the  smoothing 
by  requiring  the  node  at  -1  to  equal  the  node  at  +1  placement. 

A  modification  to  the  smoothing  process  must  be  made  if  the  sample 
contains  zero  values  of  h.  In  the  upper  part  of  Figure  1,  it  will  be 
noted  that  the  node  at  3  may  be  greater  than  or  less  than  100  percent 
probability.  If  this  node  is  greater  than  or  equal  to  100,  then  the 
sample  does  not  contain  zeroes.  The  30  classes  will  contain  all  the 
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values  in  the  sample  and  will  reach  100  percent.  The  point  where  the 
smoothing  curve  crosses  the  100  percent  line  can  be  reflected  back  to  h 
and  is  an  estimate  of  the  lowest  value  of  h. 

If  the  sample  contains  zeroes  and  the  zeroes  are  excluded  from  the 
30-class  tally,  the  sample  probabilities  will  not  reach  100  percent, 
and  node  three  will  be  less  than  100.  In  this  case,  the  zeroes  are 
considered  to  compose  a  31st  class,  and  provide  a  sample  point  on  the 
100  percent  probability  line.  This  requires  an  extra  node  placed  at  5, 
and  also  requires  an  iterative  least  squares  smoothing  to  find  where  a 
fictitious  arc  crosses  the  100  percent  line.  This  arc  is  called  ficti¬ 
tious  because  it  implies  values  of  h  less  than  zero,  and  this  is 
impossible.  The  method  has  the  advantage,  however,  of  providing  analy¬ 
tical  integrity  across  the  non-zero  and  zero  items  in  a  mixed  sample. 

Program  1,  below,  performs  smoothing  of  data  without  zeroes. 
Program  2  performs  smoothing  of  data  with  zeroes.  Note  that  no 
distinction  is  made  between  extreme-value  maxima,  extreme-value  minima, 
nor  total  samples. 

Further  work  may  be  needed  on  the  best  value  of  the  empirical 
parameter  k.  It  could  possibly  be  a  function  of  the  third  moment  of 
the  sample.  The  transform  from  h  to  v  may  also  need  some  relatively 
minor  modifications.  In  particular  the  common  point  of  the  two  ex¬ 
ponential  limbs  may  need  to  be  moved  from  v=1.5  toward  v=2.0  for  some 
highly  skewed  distributions  of  maxima.  Also,  the  calibration  point 
could  be  changed  from  h  1  = - 1 . 0 ,  v  =  3.0.  The  only  objective  of  the  h  to 
v  transform  is  to  produce  a  fairly  smooth  sample  probability  rising 
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METHOD  OF  SIMULATION 

Programs  1  and  2  analyze  one  historical  sample.  The  derived 
smoothing  curves  are  estimates  of  population  probabilities  as  provided 
by  the  sample.  We  often  need  estimates  of  possible  variability  if 
such  samples  were  to  be  repeated.  For  example,  we  might  be  interested 
in  possible  future  daily  or  monthly  rainfall  variability.  We  might  be 
interested  in  how  many  times  a  critical  value  could  be  equalled  or 
exceeded  in,  say,  100  such  future  samples.  These  values  form  confidence 
intervals  and  tolerance  limits  on  values  for  planning  and  design. 

The  smoothing  curves  derived  from  samples  by  Program  1  and  Program 
2  can  be  used  to  synthesize  possible  future  samples.  Random  numbers 
can  be  generated  in  a  computer  program.  Each  such  random  number,  in  a 
scale  from  zero  to  100,  may  be  considered  a  probability.  Each  value 
of  probability  has  an  equal  chance  of  occurrence.  These  random  proba¬ 
bilities  can  be  reflected  through  the  smoothing  curves  by  reverse 
interpolation  to  yield  values  of  v.  These  random  v-values  can  then 
be  transformed  to  random  h's.  As  an  example,  if  h  represented  the 
annual  maximum  flood,  then  50  random  h's  could  be  called  a  simulated 
50-year  record  of  annual  floods. 

Program  3,  below,  performs  simulation  of  samples  with  no  zeroes. 

The  smoothing  curve  would  have  been  derived  using  Program  1  in  the 
analysis  of  an  historical  sample. 

Program  4  performs  simulation  of  samples  with  zeroes.  Program  2 
would  have  been  used  for  analysis. 

Both  programs  allow  for  simulation  of  a  varying  number  of  simulated 
samples  of  a  designated  number  of  items  each.  Again,  an  example,  if 
h  represented  the  annual  maximum  flood,  100  sets  of  50  random  h's  could 
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be  called  100  possible  future  50-year  simulated  flood  records.  Each  of 
the  100  are  placed  in  order  of  magnitude.  For  rank  number  1  of  the  50, 
one  has  thus  a  100  item  distribution  of  the  l-in-50  year  flood.  The 
1  - in-50  sets  a  risk  level.  The  100  item  distribution  is  the  uncertainty 
associated  with  the  designated  risk. 
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Figure  1.  Geometrical  Schematic  of  h  to  v  Transform  for  Placement 
of  Nodes. 
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APPENDIX 


The  Appendix  includes  four  program  listings  and,  notes  and  sample 
outputs  of  each.l/ 


-^These  programs  are  presented  only  for  the  convenience  of  potential 
users.  While  the  programs  have  been  run  and  tested  on  various 
data  sets  the  originators  of  the  programs  assume  no  responsibility 
for  either  their  accuracy  or  adequacy.  Such  responsibility  must 
rest  solely  on  the  user.  We  stand  ready  to  assist  and  advise 
within  the  restrictions  imposed  by  our  operating  resources. 
Programs  1  and  2  are  written  in  CYBER  BASIC.  Programs  3  and  4  are 
written  in  FORTRAN  G  (Trade  name  is  included  for  the  benefit  of 
the  reader  and  does  not  imply  an  endorsement  or  preferential 
treatment  of  the  named  products). 
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PROGRAM  1 


1 1 


To  be  used  on  data  sets  with  few  zeroes.  Node  No.  3  must  be  _>  1 00 
probability.  (See  Figure  1). 


This  program  finds  four  nodes  of  the  sliding  polynomials  with  data 
organized  into  30  classes. 


Boundary  condition  is  set  so  that  the  sliding  polynomial  smoothing 
curve  is  zero  with  zero  slope  at  +  °°  of  the  data  scale. 


Use  Program  3  for  simulation  with  the  nodes  derived  from  this  fitting. 
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00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 


REM  ANALYSIS  WITH  FEW  ZEROES 
REM  READ  FROM  FILE  //I  AND  WRITE  TO  FILE  #5 
OPTION  BASE  1 
FILE  #1="PER1" 

RESTORE  in 
FILE  #5="SMDATA" 

RESTORE  #5 

DIM  H(200) ,W(4,4) ,M(35) , C ( 35 ) , P(35) ,Y(6) 

DIM  R ( 3 5 ) ,X(35) , T(35 ) , V ( 35  ) , S (4) 

DIM  U ( 35 ) 

REM  INITIALIZE  ARRAYS 
FOR  1=1  TO  30 
NODATA  00290 

READ  R ( I ) ,X(I) ,T(  I ) , V( I ) 

NEXT  I 

DATA  0.028,-0.0045,0,0,0. 104,-0.016,0,0,0.216,-0.0315,0,0 
DATA  0.352,-0.048,0,0,0.5,-0.0625,0,0,0.648,-0.072,0,0 
DATA  0.784,-0.0735,0,0,0.896,-0.064,0,0,0.972,-0.0405,0,0 
DATA  1,0,0,0,0.9765,0.0685,-0.0045,0,0.912,0. 168,-0.016,0 
DATA  0.8155,0.2895,-0.0315,0,0.696,0.424,-0.048,0,0.5625 
DATA  0.5625 , -0.0625,0,0.424,0.696,-0.072 

DATA  0,0.289  5,0.8155,-0.0735,0,0. 168,0.912,-0.064,0,0.0685 
DATA  0.9765,-0. 0405,0,0, 1 ,0 ,0 , -0 . 0405,0.9765,0. 0685 , -0 . 0045 
DATA  -0.064,0.912,0.168,-0.016,-0.0735,0.8155,0.2895,-0.0315 
DATA  -0.072,0.696,0.424,-0.048,-0.0625,0.5625,0.5625 
DATA  -0.0625,-0.048,0.424,0.696,-0.072,-0.0315,0.2895,0.8155 
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00270  DATA  -0,0735,-0.016,0.168,0.912,-0.064,-0.0045,0.0685,0.9765  27 
00280  DATA  -0.0405,0,0,1,0  28 
00290  FOR  1=1  TO  4  29 
00300  FOR  J=1  TO  4  30 
00310  READ  W(I,J)  31 
00320  NEXT  J  32 
00330  NEXT  I  33 
00340  DATA  0.133799,-0.0277062,0.00404919,-0.261391  34 
00350  DATA  -0.0277062,0.157706,0.0568727,1.86251  35 
00360  DATA  0.00404919,0.0568727,0.746445,9.06661  36 
00370  DATA  -0.261391,1.86251,9.06661,161.165  37 
00380  PRINT  #5,,  "DEFICIENT  RAIN  PERIOD  #l"  38 
00390  :MV  VM-MV  MV  MV.  MU  39 
00400  :VMVV  40 
00410  PRINT  #5,  41 
00420  PRINT  #5,  "INPUT  HYDROLOGIC  SAMPLE”  42 
00430  PRINT  #5  43 
00440  REM  INPUT  N  NUMBER  OF  DATA  44 
00450  INPUT  #1 ,N  45 
00460  PRINT  #5 ,  "NUMBER  OF  ITEMS  IN  SAMPLE. ";N  46 
00470  PRINT  #5  47 
00480  REM  PARAMETER  K=2.0  48 
00490  INPUT  V 1,K  49 
00500  REM  INPUT  DATA  INTO  H(  )  50 
00510  FOR  1=1  TO  N  51 
00520  INPUT  #  1 ,  H  ( I )  52 
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00530  PRINT  #5  USING  00400,  H(I);  53 
00540  NEXT  I  54 
00550  PRINT  U 5  55 
00560  S1=0  56 
00570  S3=0  57 
00580  REM  S 1=SUM  OF  DATA,  S3=SUM  OF  SQUARES  58 
00590  FOR  1=1  TO  N  59 
00600  S 1=S1+H ( I )  60 
00610  S 3=S 3+H ( I ) *H ( I )  61 
00620  NEXT  I  62 
00630  H1=S1/N  63 
00640  S5=SQR( (S3-S1*S1/N)/ (N-l) )  64 
00650  PRINT  #5  65 
00660  PRINT  #5,  "  SAMPLE  AVERAGE  IS";H1  66 
00670  PRINT  #5  67 
00680  PRINT  #5,  "  SAMPLE  SD  IS  ";S5  68 
00690  1=0  69 
00700  S4=K”-S5  70 
00710  V0=4-2 . 5 ‘'-'EXP (-0.9  1629-H1/S4)  71 
00720  V0=INT ( 100*V0+0 . 5 ) / 100  72 
00730  V2=V0/30  73 
00740  FOR  V 1=V2  TO  VO  STEP  V2  74 
00750  1=1+1  75 
00760  IF  V 1 > 1 . 5  THEN  00800  76 
00770  REM  M (  )=CLASS  LIMITS  77 
00780  M ( I )=LOG (VI/ 1 . 5 )*S4/ ( - 1 . 527 15 )+Hl  78 
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00790  GOTO  00810  79 
00800  M ( I )=LOG ( (4-Vl)/2.5 )*S4/0 . 9 1629+H1  80 
00810  NEXT  VI  81 
00820  FOR  L=1  TO  35  82 
00830  C(L)=0  83 
00840  NEXT  L  84 
00850  REM  COMPUTE  CLASS  FREQUENCY,  C(  )  85 
00860  FOR  1=1  TO  N  86 
00870  K=1  87 
00880  IF  H(I )<M (K)  THEN  00910  88 
00890  C(K)=C(K)+1  89 
00900  GOTO  00930  90 
00910  K=K+1  gi 
00920  GOTO  00880  92 
00930  NEXT  I  93 
00940  FOR  L=1  TO  30  94 
00950  REM  U(  )=CLASS  SAMPLE  95 
00960  U(L)=C(L)  96 
00970  NEXT  L  97 
00980  A1=0  98 
00990  REM  COMPUTE  CLASS  PROBABILITY,  P(  )  99 
01000  FOR  1=1  TO  30  100 
01010  A1=A1+C ( I )  10i 
01020  P ( I ) =A 1 / N  102 
01030  NEXT  I  103 
01040  FOR  1=1  TO  4  104 
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01050 

s (I )=0 

105 

01060 

NEXT  I 

106 

01070 

FOR  1=1  TO  30 

107 

01080 

S(1)=S(1)+R(I)*P(I) 

108 

01090 

S(2)=S(2)+X(I)*P(I) 

109 

01100 

S(3)=S(3)+T(I)*P(I) 

110 

onio 

S(4)=S(4)+V(I)*P(I) 

111 

01120 

NEXT  I 

112 

01130 

FOR  1=1  TO  4 

113 

01140 

Y ( I ) =0 

114 

01150 

FOR  J=1  TO  4 

115 

01160 

Y(I)=Y(I)+W(I,J)*S(J) 

116 

01170 

NEXT  J 

117 

01180 

NEXT  I 

118 

01190 

PRINT  #5 

119 

01200 

PRINT  #5 ,  "  SLIDING  POLYNOMIAL  ORDINATES" 

120 

01210 

FOR  L=1  TO  4 

121 

01220 

PRINT  it 5  ,  Y  (L) 

122 

01230 

NEXT  L 

123 

01240 

FOR  L=1  TO  4 

124 

01250 

K=5  -L 

125 

01260 

Y (K+2 )=Y (K) 

126 

01270 

PRINT  #5 

127 

01280 

NEXT  L 

128 

01290 

Y ( 2 )=0 

129 

01300 

Y  ( 1 )  =  Y  ( 3 ) 

130 
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01310  REM  COMPUTE  SMOOTH  PROBABILITY,  C(  )  131 
01320  FOR  K=2  TO  4  132 
01330  Z=-0.5  133 
01340  A=(9* (Y(K)+Y(K+1 ) ) -Y(K- 1 ) -Y (K+2) ) / 16  134 
01350  B=(ll*(Y(K+l)-Y(K))+Y(K-l)-Y(K+2))/8  135 
01360  Q=(Y(K-1)-Y(K)-Y(K+1)+Y (K+2 ) ) /4  136 
01370  D=(3*(Y(K)-Y(K+1))-Y(K-1)+Y (K+2) ) /2  137 
01380  FOR  J=1  TO  10  138 
01390  Z=Z+0.1  139 
01400  I 1=(K-2)*10+J  140 
01410  C(I1)=((D*Z+Q)*Z+B)*Z+A  141 
01420  NEXT  J  142 
01430  NEXT  K  143 
01440  PRINT  #5,"  CLASS  CLASS  SAMPLE  SMOOTH"  144 
01450  PRINT  #5,"  LIMITS  SAMPLE  PROBABILITY  PROBABILITY"  145 
01460  PRINT  #5  146 

01470  PRINT  # 5,"  . . . ---- . .  -------- . "  147 

01480  FOR  1=1  TO  30  148 
01490  PRINT  # 5  USING  00390 , I , M ( I ) , U( I ) , P ( I ) , C ( I )  149 
01500  NEXT  I  150 
01510  END  151 
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Notes  for  Program  1 


Line  # 

3 

16-28 

34  -  38 

45  -  54 
56  -  68 
69  -  73 
74  -  81 
86  -  93 
100-103 
107-112 
113-118 

121-130 

132-143 


Comment 


Begins  arrays  at  one. 

Sliding  Polynomial  Coefficients  for  data  set 
organized  into  30  classes. 

Inverse  of  Sums-of-Products  matrix  of  Sliding 
Polynomial  Coefficients.  This  matrix  is  fixed 
for  30  classes  and  4  nodes. 

Input  Sample. 

Compute  average  and  standard  deviation. 

Set  width  of  30  uniform  classes  in  v-scale. 

Compute  class  boundaries  in  original  data  scale. 

Tally  the  sample  into  the  30  classes. 

Accumulate  across  classes  and  take  ratios. 

Compute  the  z  XY  vector  of  least  squares. 

Multiply  Z  XY  vector  by  inverse  matrix  to  get  4 
nodes . 

Re-position  the  nodes  and  put  in  the  boundary 
nodes . 

Lay  in  the  ensemble  of  3  Sliding  Polynomial  arcs 
across  4  nodes. 
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Sample  Output  For  Program  1 


Deficient  Rain  Pericd  #1 
Input  Hydrclogic  Sample 


Number  of  Items  in  Sample  116 


21 

25 

5 

7 

6 

44 

11 

10 

6 

36 

19 

8 

1 

9 

8 

12 

20 

3 

3 

12 

17 

39 

5 

2 

13 

11 

15 

6 

32 

30 

28 

2 

1 

13 

21 

27 

5 

12 

42 

15 

7 

31 

29 

11 

36 

22 

5 

18 

25 

31 

3 

4 

1 

2 

10 

8 

14 

1 

9 

18 

5 

1 

2 

6 

28 

11 

9 

2 

5 

28 

3 

1 

21 

12 

7 

13 

9 

2 

24 

1 

11 

23 

27 

13 

43 

25 

9 

23 

9 

7 

5 

36 

7 

29 

8 

1 

3 

7 

16 

6 

12 

6 

28 

8 

28 

35 

30 

12 

19 

36 

2 

3 

32 

24 

5 

50 

Sample  Average  is  14.7845 

Sampl e  SD  is  11 .8958 

Sliding  Polynomial  Ordinates 

0.259249 

0.442444 

1 .00052 

1 .86506 


Cl  ass 

Cl  ass 

Sampl  e 

Smooth 

Limits 

Sampl e 

Probabi 1 i ty 

Probabi 1 

1 

59.263 

0 

0.00000 

0.0053 

2 

48.464 

1 

0.00862 

0.0199 

3 

42.147 

2 

0.02586 

0.0421 

4 

37.666 

2 

0.04310 

0.0700 

5 

34.189 

5 

0.08621 

0.1020 

6 

31 .349 

2 

0.10345 

0.1361 

7 

28.947 

6 

0.15517 

0.1707 

8 

26.867 

7 

0.21552 

0.2040 

9 

25.032 

0 

0.21552 

0.2341 

10 

23.391 

5 

0.25862 

0.2592 

11 

21.906 

3 

0.28448 

0.2790 

12 

20.550 

3 

0.31034 

0.2948 

13 

19.303 

1 

0.31897 

0.3080 

14 

18.149 

2 

0.33621 

0.3200 

15 

17.074 

2 

0.35345 

0.3322 

16 

16.068 

1 

0.36207 

0.3458 

17 

15.124 

1 

0.37069 

0.3623 

18 

14.217 

2 

0.38793 

0.3830 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 


Sample  Output  For  Program  1  (continued) 


Class 

Cl  ass 

Sampl  e 

Smooth 

Limi ts 

Samp! e 

Probabi  1  i  t.y 

Probabi 1 i 

13.284 

1 

0.39655 

0.4093 

12.317 

4 

0.43103 

0.4424 

11  .311 

6 

0.48276 

0.4817 

10.266 

5 

0.52586 

0.5252 

9.176 

2 

0.54310 

0.5727 

8.039 

6 

0.59483 

0.6240 

6.849 

11 

0.68966 

0.6789 

5.602 

6 

0.74138 

0.7372 

4.293 

8 

0.81034 

0.7988 

2.914 

7 

0.87069 

0.8633 

1  .457 

7 

0.93103 

0.9306 

-0.086 

8 

1 .00000 

1 .0005 

Figure  2  is  a  plot  of  this  sample  output 


PROBABILITY  OF  N  >  SCALE  N 


NUMBER  OF  DAYS  -  CLASS  LIMITS 


Figure  2.  Distribution  of  Number  of  Sequential  Days  with  Rain 
Less  Than  1 .0  cm. 
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PROGRAM  2 


To  be  used  on  data  sets  containing  zeroes. 


This  program  finds  five  nodes  of  the  sliding  polynomials.  Non-zero 
data  are  organized  into  30  classes.  Zero  data  values  are  put  into  a 
31st  class.  This  program  has  performed  satisfactorily  for  the  authors 
with  data  set  containing  as  much  as  30  percent  zeroes. 


Left  boundary  condition  is  set  so  that  the  sliding  polynomial  smoothing 
curve  is  zero  with  zero  slope  at  +  °°  of  the  data  scale.  Right  boundary 
is  free.  Difference  between  100  and  node  number  4  is  estimate  of  popu¬ 
lation  proportion  of  zeroes.  A  fifth  node  is  necessary  and  is  assumed 
to  lie  on  a  parabola  through  nodes  2,  3,  and  4. 


Use  Program  4  for  simulation  with  the  nodes  derived  from  this  fitting. 


23 


00010  REM  ANALYSIS  WITH  ZEROES  1 
00020  REM  READ  FROM  FILE  #1  AND  WRITE  INTO  FILE  #5  2 
00030  OPTION  BASE  1  3 
00040  FILE  #1=MJUL"  4 
00050  RESTORE  #1  5 
00060  FILE  #5="SMDATA"  *  6 
00070  RESTORE  #5  7 
00080  DIM  H(2000) , M(35 ) ,0(35) ,  P (4 , 3 1 )  8 
00090  DIM  E (4 , 4) , G (4 , 4) , Q (4 , 4) ,Y(7)  9 
00100  DIM  A$(80) ,R(4) ,D(35) ,T(40)  10 
00110  FOR  1=1  TO  30  11 
00120  FOR  J=1  TO  4  12 
00130  NODATA  00300  13 
00140  READ  P(J,I)  14 
00150  NEXT  J  15 
00160  NEXT  I  16 
00170  DATA  0.028,-0.0045,0,0,0.104,-0.016,0,0,0.216,-0.0315,0,0  17 
00180  DATA  0.352,-0.048,0,0,0.5,-0.0625,0,0,0.648,-0.072,0,0  18 
00190  DATA  0.784,-0.0735,0,0,0.896,-0.064,0,0,0.972,-0.0405,0,0  19 
00200  DATA  1,0,0,0,0.9765,0.0685,-0.0045,0,0.912,0.168,-0.016,0  20 
00210  DATA  0.8155,0.2895,-0.0315,0,0.696,0.424,-0.048,0,0.5625  21 
00220  DATA  0.5625,-0.0625,0,0.424,0.696,-0.072  22 
00230  DATA  0,0.2895,0.8155,-0.0735,0,0.168,0.912,-0.064,0,0.0685  23 
00240  DATA  0 . 9765 , -0 . 0405 , 0 , 0 , 1 , 0 , 0 , -0 . 0405 , 0 . 9 765 , 0 . 0o85 , -0 . 0045  24 
00250  DATA  -0.064,0.912,0.168,-0.016,-0.0735,0.8155,0.2895,-0.0315  25 
00260  DATA  -0.072,0.696,0.424,-0.048,-0.0625,0.5625,0.5625  26 
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00270  DATA  -0.0625,-0.048,0.424,0.696,-0.072,-0.0315,0.2895,0.8155  27 

00280  DATA  -0.0735,-0.016,0.168,0.912,-0.064,-0.0045,0.0685,0.9765  28 

00290  DATA  -0.0405,0,0,1,0  29 

00300  FOR  1=1  TO  4  30 

00310  FOR  J=1  TO  4  31 

00320  E ( I , J)=0  32 

00330  NEXT  J  33 

00340  NEXT  I  34 

00350  :###  ww.vnu  #####  Ur^n.nm  35 

00360  '.mm  mm  36 

00370  -.pmimm  37 

00380  REM  PRINT  OPTION  IF  A# 1 , ALL  IS  PRINTED  38 

00390  A=1  39 

00400  PRINT  #5,  --------------JUL . . . "  40 

00410  PRINT  #5  41 

00420  PRINT  ft 5,  "  — . - . "  42 

00430  PRINT  #5  43 

00440  PRINT  #5,  ---------SAMPLE  ITEMS - - "  44 

00450  PRINT  ft 5  45 

00460  N=1  46 

00470  IF  END  ftl  THEN  00520  47 

00480  INPUT  #1 ,H(N)  48 

00490  PRINT  ft 5  USING  00370, H(N);  49 

00500  N=N+1  50 

00510  GO  TO  00470  51 

00520  N=N - 1 


52 
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00530  PRINT  #5  53 
00540  K=2 . 0  54 
00550  S3=0  55 
00560  S1=0  56 
00570  FOR  1=1  TO  N  57 
00580  S 1=S 1+H ( I )  58 
00590  S3=S3+H ( I )*H ( I )  59 
00600  NEXT  I  60 
00610  H1=S1/N  61 
00620  S5=SQR ( (S3-S1*S 1/N) / ( N -  1 ) )  62 
00630  PRINT  #5  63 
00640  PRINT  #5,  ”  SAMPLE  AVERAGE  IS  ";H1  64 
00650  PRINT  #5,  "  SAMPLE  STANDARD  DEVIATION  IS  ";S5  65 
00660  PRINT  #5  66 
00670  PRINT  # 5,  "---CLASS  LIMITS -- -CLASSIFIED  SAMPLE - -CLASS  PROBABILITI  67 
00680  1=0  68 
00690  V0=4-2.5*EXP(-.91629*H1/K/S5)  ’  69 
00700  V0=INT ( 1000*V0+0 . 5 ) / 1000  70 
00710  V2=V0/30  71 
00720  FOR  V1=V2  TO  VO  STEP  V2  72 
00730  1=1+1  73 
00740  IF  Vl>=  1.5  THEN  00770  74 
00750  M(I)=K*S5*L0G(V1/1. 5)/ (- 1.52715 )+Hl  75 
00760  GO  TO  00780  76 
00770  M ( I )=K*S5*LOG ((4-Vl)/2.5)/0. 91629 +H1  77 


00780  NEXT  VI 


78 
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00790  FOR  L=1  TO  31 
00800  C (L)=0 
00810  NEXT  L 
00820  1=1 

00830  IF  H ( I )=0  THEN  00870 
00840  K=1 

00850  IF  H(I )>M (K)  THEN  00910 

00860  GO  TO  00890 

00870  C (31 )=C (31 )+l 

00880  GO  TO  00920 

00890  K=K+1 

00900  GO  TO  00850 

00910  C (K)=C (K)+l 

00920  IF  I=N  THEN  00950 

00930  1=1+1 

00940  GO  TO  00830 

00950  FOR  1=1  TO  31 

00960  D(I)=C(I) 

00970  NEXT  I 
00980  A1=0 
00990  FOR  1=1  TO  31 
01000  A1=A 1+C ( I ) 

01010  C(I)=A1/N 

01020  PRINT  #5  USING  00350 , 1 , M ( I ) , D ( I )  , C 
01030  NEXT  I 
01040  1=1 


79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 

(I)  102 

103 

104 
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01050  K=1 

105 

01060  L=1 

106 

01070  E (I ,K)=E(I,K)+P(I,L)*PCK,L) 

107 

01080  IF  L=30  THEN  OHIO 

108 

01090  L=L+1 

109 

01100  GOTO  01070 

110 

OHIO  IF  K=4  THEN  01140 

111 

01120  K=K+1 

112 

01130  GOTO  01060 

113 

01140  IF  1=4  THEN  01170 

114 

01150  1=1+1 

115 

01160  GOTO  01050 

116 

01170  1=1 

117 

01180  FOR  K=1  TO  4 

118 

01190  E(K,I)=E(I,K) 

119 

01200  NEXT  K 

120 

01210  PRINT  #5 

121 

01220  Z=0 . 25 

122 

01230  PRINT  #5,  "--INITIAL  Z  VALUE =  ";Z 

123 

01240  S2=Z 

124 

01250C1=(((-8*Z+4)*Z+2)*Z-1)/16 

125 

01260  C2=( ( (24*Z-4)*Z-22)*Z+9)/16 

126 

01270  C3= ( ( ( - 24*Z - 4 ) *z+2 2 ) *Z+9 ) / 1 6 

127 

01280  C4=(((8*Z+4)*Z-2)*Z-1)/16 

128 

01290  IF  A=1  THEN  01330 

129 

01300  PRINT  #5 

130 

28 

01310  PRINT  #5,  "CLASS  31  COEFS .  " ; C 1 ; C2 ; C3 ; C4  131 

01320  PRINT  #5  132 

01330  P ( 1 , 3 1 )=0  133 

01340  P (2 , 3 1 )=C 1+C4  134 

01350  P(3,31)=C2-3*C4  135 

01360  P  (4 , 3 1  )=C3+3'VC4  136 

01370  G ( 2 , 2 ) =E ( 2 , 2 ) +P ( 2 , 3 1 ) *P ( 2 , 3 1 )  137 

01380  G ( 2 , 3 ) =E ( 2 , 3 ) +P ( 2 , 3 1 ) *P ( 3 , 3 1 )  138 

01390  G(3,2)=G(2,3)  139 

01400  G(2,4)=E(2,4)+P(2,31)*P(4,31)  140 

01410  G(4,2)=G(2,4)  141 

01420  G(3,3)=E(3,3)+P(3,31)*P(3,31)  142 

01430  G(3,4)=E(3,4)+P(3,31)*P(4,31)  143 

01440  G(4,3)=G(3,4)  144 

01450  G(4,4)=E(4,4)+P(4,31)*P(4,31)  145 

01460  G ( 1 , 1 ) =E (1,1)  146 

01470  G ( 1 , 2 ) =E (1,2)  147 

01480  G ( 1 , 3 )=E (1,3)  148 

01490  G(1,4)=E(1,4)  149 

01500  G(2, 1 ) =E (2,1)  150 

01510  G ( 3 , 1 ) =E (3,1)  151 

01520  G(4, 1)=E(4, 1)  152 

01530  IF  A=1  THEN  01620  153 

01540  PRINT  #5 ,  " - G(I,J)  MATRIX----"  154 

01550  FOR  1=1  TO  4  155 

01560  PRINT  #5,  I  156 
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01570  FOR  J=1  TO  4 

157 

01580  PRINT  #5 ,  G(I,J); 

158 

01590  NEXT  J 

159 

01600  PRINT  #5 

160 

01610  NEXT  I 

161 

01620  FOR  1=1  TO  4 

162 

01630  FOR  J=1  TO  4 

163 

01640  Q ( I , J)=0 

164 

01650  NEXT  J 

165 

01660  NEXT  I 

166 

01670  FOR  1=1  TO  4 

167 

01680  Q ( I , I )=1 

168 

01690  NEXT  I 

169 

01700  J=1 

170 

01710  I=J 

171 

01720  IF  G ( I , J)<>0  THEN  01780 

172 

OT 7  30  IF  1=4  THEN  01760 

173 

01740  1=1+1 

174 

01750  GOTO  01720 

175 

01760  PRINT  #5,  "SINGULAR  MATRIX" 

176 

01770  GOTO  02420 

177 

01780  GOTO  01790 

178 

01790  K=1 

179 

01800  S=G ( J , K) 

180 

01810  G ( J , K ) =G ( I , K ) 

181 

01820  G  ( I , K)=S 

182 

30 

01830  S=Q ( J , K) 

183 

01840  Q(J,K)=Q(I,K) 

184 

01850  Q(I ,K)=S 

185 

01860  IF  K=4  THEN  01890 

186 

01870  K=K+1 

187 

01880  GOTO  01800 

188 

01890  T=1/G(J,J) 

189 

01900  K=1 

190 

01910  G(J,K)=T*G(J,K) 

191 

01920  Q ( J , K)=T*Q ( J  ,  K) 

192 

01930  IF  K=4  THEN  01960 

193 

01940  K=K+1 

194 

01950  GOTO  01910 

195 

01960  L=1 

196 

01970  IF  L=J  THEN  02050 

197 

01980  T=-G(L,J) 

198 

01990  K=1 

199 

02000  G(L,K)=G(L,K)+T*G(J,K) 

200 

02010  Q(L,K)=Q(L,K) +T*Q ( J , K ) 

201 

02020  IF  K=4  THEN  02050 

202 

02030  K=K+ 1 

203 

02040  GOTO  02000 

204 

02050  IF  L=4  THEN  02080 

205 

02060  L=L+1 

206 

02070  GOTO  01970 

207 

02080  IF  J=4  THEN  02110 

208 

31 

02090  J=J+1 

209 

02100  GOTO  01710 

210 

02110  IF  A=1  THEN  02210 

211 

02120  PRINT  #5 

212 

02130  PRINT  #5,  - -INVERSE  MATRIX--- 

ft 

213 

02140  FOR  1=1  TO  4 

214 

02150  PRINT  #5,  I 

215 

02160  FOR  J=1  TO  4 

216 

02170  PRINT  #5,  Q(I,J); 

217 

02180  NEXT  J 

218 

02190  PRINT  #5 

219 

02200  NEXT  I 

220 

02210  R ( 1 )=0 

221 

02220  R(2)=0 

222 

02230  R (3 )=0 

223 

02240  R (4 )=0 

224 

02250  1=1 

225 

02260  R ( 1 )  =R  ( 1 )  +P ( 1 , 1 ) ••C  ( I ) 

226 

02270  R(2)=R(2)+P(2,I)*C(I) 

o  9  7 

02280  R ( 3 ) =R ( 3 ) +P ( 3 , I ) *C ( I ) 

228 

02290  R ( 4 ) =R ( 4 ) +P ( 4 , I ) *C ( I ) 

229 

02300  IF  1=31  THEN  02330 

230 

02310  1=1+1 

231 

02320  GOTO  02260 

232 

02330  1=1 

233 

02340  Y( I )=0 

234 

32 

02350  J=1  235 
02360  Y ( I ) =Y ( I ) +Q ( I , J ) *R ( J )  236 
02j70  IF  J=4  THEN  02400  737 
02380  J=J+1  938 
02390  GOTO  02360  239 
02400  IF  1=4  THEN  02430  240 
02410  1=1+1  241 
02420  GOTO  02340  242 
02430  Y(5)=Y(2) -3*Y (3 )+3*Y (4)  243 
02440  IF  A=1  THEN  02480  244 
02450  PRINT  #5,  ----------SOLUTION  NODES----------  245 

02460  PRINT  #5,  Y ( 1 ) ;Y(2) ; Y ( 3 ) ;Y(4) ;Y(5)  246 
02470  PRINT  #5  247 
02480  Z l=-0 . 5  248 
02490  Z2=0 . 5  249 
02500  Z=(Zl+Z2)/2  250 
02510  C 1=  ( ( ( - 8* Z+4 )'•'■  Z+2 )* Z - 1 )  /  16  251 
02520  C2=(((24*Z-4)*Z-22)*Z+9)/16  252 
02530  C3=(((-24*Z-4)*Z+22)*Z+9)/16  253 
02540  C4=(  ( ( 8*Z+4) *Z - 2 )'“Z -  1 )  / 16  254 
02550  P1=C1*Y(2)+C2*Y(3)+C3*Y(4)+C4*Y(5)  255 
02560  IF  A=1  THEN  02580  256 

02570  PRINT  # 5,  " - PI  NOW  IS  ";P1  257 

02580  IF  ABS (PI- 1)<0 . 001  THEN  02640  258 
02590  IF  P1>1  THEN  02620  259 
02600  Z1=Z  260 
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02610  GOTO  02500  261 
02620  Z2=Z  262 
02630  GOTO  02500  263 
02640  IF  A=1  THEN  02670  264 
02650  PRINT  #5  265 
02660  PRINT  #5,  "VALUE  OF  100%  INTERCEPT  IS  ";3.5+Z  266 
02670  IF  ABS (S2-Z) <0.001  THEN* 02700  267 
02680  S2=Z  268 
02690  GOTO  01250  269 
02700  PRINT  #5  270 
02710  PRINT  #5 ,  "  CONVERGENCE  STOP  "  271 
02720  PRINT  #5  272 
02730  PRINT  #5 , ------SOLUTION  NODES----"  273 

02740  PRINT  #5,  Y ( 1 ) ; Y ( 2 ) ; Y ( 3 ) ; Y (4 ) ; Y (5 )  274 
02750  REM  RELOCATE  ORDINATES  275 
02760  Y(7)=Y(5)  276 
02770  Y(6)=Y (4)  277 
02780  Y(5)=Y(3)  278 
02790  Y(4)=Y(2)  279 
02800  Y(3)=Y(1)  280 
02810  Y(2 )=0  281 
02820  Y ( 1 ) =Y ( 3 )  282 
02830  FOR  K=2  TO  5  283 
02840  Z=-0 . 5  284 
02850  A= ( 9* ( Y (K ) +Y (K+ 1 ) ) - Y (K - 1 ) - Y (K+2 ) ) / 1 6  285 
02860  B=( 1 1* ( Y (K+l ) -Y (K) )+Y(K- 1 ) -Y (K+2) ) /8  286 
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02870  C=(Y(K-l)-Y(K)-Y(K+l)+Y(K+2))/4  287 
02880  D=(3* (Y(K) -Y (K+l) ) -Y(K- 1 )+Y (K+2) ) /2  288 
02890  FOR  J=1  TO  10  289 
02900  Z=Z+0 . 1  °90 
02910  I 1=(K-2)*10+J  291 
02920  T ( 1 1 )  =  ( ( D* Z +C ) * Z+B ) * Z + A  292 
02930  NEXT  J  293 
02940  NEXT  K  294 
02950  PRINT  ,45  295 
02960  PRINT  #5,  "  INTERPOLATED  VALUES  "  296 
02970  PRINT  #5  297 
02980  FOR  1=1  TO  40  298 
02990  PRINT  45,  I ;T( I ) ,  299 
03000  NEXT  I  300 
03010  PRINT  45  30i 
03020  STOP  302 
03030  END  303 
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Notes  For  Program  2 

Line  # 

Comment 

3 

Begins  arrays  at  one. 

17  -  29 

Sliding  Polynomial  Coefficients  for  data  set 
organized  into  30  classes. 

47  -  52 

Input  sample. 

55  -  65 

Compute  average  and  standard  deviation. 

68  -  71 

Set  width  of  30  uniform  classes  in  v-scale. 

72  -  78 

Compute  class  boundaries  in  original  data  scale. 

82  -  94 

Tally  the  non-zero  values  of  the  sample  into  30 
classes.  Zeroes  tally  into  class  31. 

98-103 

Accumulate  across  classes  and  take  ratios. 

104  -  120 

Compute  the  "Sums-of-Products"  matrix  for  30  classes 
of  non-zeroes. 

122  -  128 

Initialize  z  of  31st  class  with  trial  value  of  0.25 
and  compute  Sliding  Polynomial  Coefficients  of  31st 
class. 

133  -  152 

Augment  the  fixed  "Sums-of-Products"  matrix  of  30 
classes  with  Sums-of-Products  of  trial  Coefficients 
of  31st  class.  Fifth  node  is  assumed  on  extension 
of  parabola  through  nodes  2,  3,  &  4. 

162  -  210 

Invert  the  Sums-of-Products  matrix. 

221  -  232 

Compute  the  E  XY  vector  of  least  squares. 

233  -  242 

Multiply  E  XY  vector  by  inverse  matrix  to  get  4  nodes 

243 

Extend  parabola  to  fifth  node. 

248  -  255 

258  -  263 

Find  new  value  of  z  where  smoothing  curve  crosses 

P  =  1 00  line.  Use  reverse  interpolation  by  interval¬ 
halving  method. 

267  -  269 

If  new  z  different  from  trial  z  of  100  percent 
crossing,  repeat  new  trial  solution. 

276  -  282 

After  iteration  to  solution,  Re-position  the  nodes 
and  put  in  the  boundary  nodes. 

283  -  294 

Lay  in  the  ensemble  of  4  sliding  polynomial  arcs 
across  5  nodes. 
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Sample  Output  For  Program  2 


July 


Sample  Items 


.1000 

2.1300 

2.6400 

3.4500 

.6900 

.9900 

.0300 

.0000 

1.6300 

.7600 

1 .0400 

.0000 

.0000 

7.8000 

2.3600 

1 .3000 

.6100 

3.5100 

4.1700 

.9700 

.3600 

3.9100 

.4100 

8.9200 

4.0100 

.0000 

1.5500 

7.2400 

.0000 

.1500 

.1500 

.0000 

.6400 

.0300 

3.4000 

.9400 

.5800 

8.3100 

2.0800 

8.1800 

1 .9600 

1 .0200 

1 .0200 

5.8900 

1 .4200 

3.3500 

4.9300 

1  .1700 

1 .4200 

.0000 

6.3800 

.4800 

.2300 

11 .2500 

.3300 

2.6200 

1  .2400 

2.3100 

11 .8400 

.4300 

.9900 

5.0800 

.0000 

5.1100 

1  .6500 

9.4200 

5.6900 

.4600 

.8600 

1 .8000 

.8100 

1 .3200 

1  .4200 

.2500 

2.7200 

4.1400 

7.1400 

8.0000 

24.9400 

.4100 

3.5600 

2.0600 

.9900 

4.2900 

1 .6800 

.0000 

2.2100 

.0000 

7.7200 

2.2900 

6.2500 

2.1100 

8.8400 

2.1300 

.0800 

.5300 

1  .1200 

.0500 

.5100 

.1000 

1 .8300 

.5300 

2.4100 

5.3100 

6.5800 

5.8900 

3.0500 

1 .5500 

1.9300 

.0500 

.0800 

2.5400 

.0500 

1 .9600 

.9400 

.8100 

.5800 

.0800 

.3000 

15.3400 

1 .5000 

3.7600 

2.9700 

5.3800 

6.8300 

.0000 

.1800 

1 .7500 

.9700 

2.1100 

.5600 

3.1500 

.0000 

.9100 

4.9500 

8.5300 

1 .3000 

2.1800 

3.2500 

.9900 

.0000 

.0000 

.4300 

.4300 

Sample  Average 

is  2.61806 

Sample  Standard 

Deviation  is 

3.40995 

Cl  ass 

Classified 

Cl  ass 

Limits 

Sampl  e 

Probabi  1  i  ties 

1 

16.01421 

1 

.00694 

2 

12.91877 

1 

.01389 

3 

11  .10806 

2 

.02778 

4 

9.82334 

0 

.02778 

5 

8.82683 

3 

.04861 

6 

8.01262 

3 

.06944 

7 

7.32422 

3 

.09028 

8 

6.72790 

3 

.11111 

9 

6.20191 

3 

.13194 

10 

5.73139 

2 

.14583 

11 

5.30576 

3 

.1  6667 

12 

4.91719 

4 

.19444 

13 

4.55973 

0 

.19444 

14 

4.22878 

1 

.20139 

15 

3.92068 

3 

.22222 
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Sample  Output  For  Program  2  (continued) 


Class 

Classified 

Cl  ass 

Limi ts 

Sampl  e 

Probabi  1  i  ties 

16 

3.63246 

2 

.23611 

17 

3.36173 

4 

.26389 

18 

3.10647 

3 

.28472 

19 

2.86502 

2 

.29861 

20 

2.63595 

2 

.31250 

21 

2.41066 

2 

.32639 

22 

2.17839 

6 

.36806 

23 

1.93865 

8 

.42361 

24 

1 .69092 

4 

.45139 

25 

1 .43466 

6 

.49306 

26 

1 .16926 

8 

.54861 

27 

.89405 

13 

.63889 

28 

.60827 

7 

.68750 

29 

.31108 

15 

.79167 

30 

.00152 

16 

.90278 

31 

.00000 

14 

1 .00000 

Initial  Z  Value  =  .25 

Convergence  Stop 

Solution  Nodes 
.154626  .322528 

.899813 

2.81714 

6. 

07451 

Interpolated  Values 

1  2.87814E-3  2 

1 . 09206E-2 

3 

2 . 32395E-2 

4 

3 . 89469E-2 

5 

5.71 549E-2 

6  7 . 69755E-2  7 

9 . 75208E-2 

8 

.117903 

9 

.137234 

10 

.1  54626 

11  .169036  12 

.180806 

13 

.191125 

14 

.20118 

15 

.212161 

16  .225254  17 

.24165 

18 

.262535 

19 

.289098 

20 

.322528 

21  .357647  22 

.390344 

23 

.423413 

24 

.459645 

25 

.501832 

26  .552766  27 

.61  5239 

28 

.692043 

29 

.785971 

30 

.899813 

31  1.03124  32 

1.17608 

33 

1.33431 

34 

1 .50594 

35 

1 .69097 

36  1.8894  37 

2.10124 

38 

2.32647 

39 

2.5651 1 

40 

2.81714 

Note:  Figure  3  is  a  plot  of  this  sample  output. 
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RAINFALL  IN  cm  -  CLASS  LIMITS 


0 


38 


igure  3.  Distribution  of  Seven-Day  Rainfall. 
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PROGRAM  3 


Simulation  of  a  number  (NS)  of  synthetic  samples  of  a  number  (NI)  of 
items  each. 


Sample  contains  few  zeroes. 


Nodes  were  derived  using  Program  1. 
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C  SIMULATION  WITH  FEW  ZEROES  --  USE  WITH  PROGRAM  NO .  1  1 

DIMENSION  H(100) , C (3 , 4) , P (6 ) , SORT (100,20) , AT ( 80 )  2 

READ (5 , 8080 )  (AT(I) ,1=1,80)  3 

8080  FORMAT (80A1)  4 

WRITE (6 , 8081)  (AT(I ) , 1=1 , 80)  5 

8081  FORMAT ( '  ' ,80A1)  6 

READ (5 , 1000 )  NR,N1 ,N2 ,NS ,NI  7 

1000  FORMAT (3 15 ,2110)  8 

WRITE (6,50)  NR,N1 ,N2 ,NS ,NI  9 

50  FORMAT ( '  ’,518)  10 

READ (5 ,1001)  (P ( I ) ,1=3,6)  11 

1001  FORMAT (4F10. 6)  12 

WRITE(6 ,8181)  (P ( I ) ,1=3,6)  13 

8181  FORMAT ( '  ' ,4F10.6)  14 

READ(5 , 1001 )  HB ,SD,Q  15 

WRITE (6 ,8181)  HB,SD,Q  16 

P ( 1 )=P (3)  17 

P(2)=0 . 0  18 

CALL  RAND (NR, N1 ,N2 ,R, 1)  19 

DO  1002  1=1,3  20 

C(I,l)=(9*(P(I+l)+P(I+2))-P(I)-P(I+3))/16  21 

C(I,2)=(ll*(P(I+2) -P(I+l))+P(I)-P(I+3))/8  22 

C(1 ,3)  =  (P(I)-P(I  +  1) -P(I+2)+P(I  +  3) )/4  23 

1002  C(I,4)=(3*(P(I+l)-P(I+2))-P(I)+P(I+3))/2  24 

VS=(4  .  -  2 . 5*EXP(-  .  9  1629V'HB/Q/SD)  )/3 . 0  25 


41 


DO  1003  M=1 , NS 

26 

DO  1004  1=1, NI 

27 

CALL  RAND ( NR , N 1 , N2 , R , 2 ) 

28 

IF  (R .  LT  .  P  ( 3 ) )  GO  TO  1005 

29 

IF  (R.LT.P(4))  GO  TO  1006 

30 

GO  TO  1007 

31 

1005 

o 

o 

II 

CO 

32 

BR=VS 

33 

I  S=  1 

34 

GO  TO  1008 

35 

1006 

BL=VS 

36 

BR=2 . 0* VS 

37 

IS=2 

38 

GO  TO  1008 

39 

1007 

BL=2 . 0*VS 

40 

BR=3 . 0*VS 

41 

IS=3 

42 

1008 

A=(»BL)/(BR-BL)-0.5 

43 

B=l/ (BR-BL) 

44 

1100 

Z=A+B* (BL+BR) /2 . 0 

45 

PC=(  (C  ( IS  ,  4)*Z+C  (IS  ,  3 )  )*Z+C  ( IS  ,  2 )  )*Z+C  ( IS  ,  1 ) 

46 

IF(R.LT.PC)  GO  TO  1009 

47 

BL=(Z-A)/B 

48 

GO  TO  1010 

49 

1009 

BR=(Z-A)/B 

50 

1010 

IF(ABS(BR-BL) .LT. 0.001)  GO  TO  1011 

51 

42 


GO  TO  1100  52 

1011  V= (BL+BR) /2  53 

IF (V . LT . 1.5)  GO  TO  1012  54 

H ( I ) =HB +Q*S  D* ALOG ((4.0-V)/2.5)/0.91629  55 

GO  TO  1004  56 

1012  H(I)=HB-Q*SD*ALOG(V/ 1.5)/ 1.52715  57 

1004  CONTINUE  58 

DO  10  J=1 , NI  59 

X6=H( J)  60 

DO  20  Jl=l ,NI  61 

IF(H(J1) .GE.X6)  GO  TO  20  62 

H(J)=H(J1)  63 

H(J1)=X6  64 

X6=H( J)  65 

20  CONTINUE  66 

10  CONTINUE  67 

WRITE (6,1017)  M  68 

1017  FORMAT (1 OX, 'SAMPLE  NO. ’,14/)  69 

WRITE(6, 1018)  (K,H(K),K=1,NI)  70 

1018  FORMAT ('  ' , 8 ( 14 , F8 . 3 ) )  71 

DO  2000  K=1 ,NI  72 

SORT(M , K)=H(K)  73 

2000  CONTINUE  74 

1003  CONTINUE  75 

DO  2010  K=1 , NI  76 
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DO  2005  1=1, NS  77 

CHECK=SORT ( I , K)  78 

DO  2004  J=1 ,NS  79 

IF(SORT(J,K) .GE. CHECK)  GO  TO  2004  80 

SORT (I ,K)=SORT( J,K)  81 

SORT(J,K)=CHECK  82 

CHECK=SORT(I,K)  83 

2004  CONTINUE  84 

2005  CONTINUE  85 

2010  CONTINUE  86 

WRITE (6 , 2223)  (L,L=1,NI)  87 

2223  FORMAT ('  ’ ,//llX, ’  SORTED  DATA  ' , / , 3X , 20 ( 12 , 4X) )  88 

DO  2001  M=1 , NS  89 

WRITE ( 6 , 2222 )  (SORT(M,K) ,K=1,NI)  90 

2222  FORMAT ( '  ' ,20F6.2)  91 

2001  CONTINUE  92 

70  STOP  93 

END  94 

SUBROUTINE  RAND (NR , N1 , N2 , DRAW , I ENT)  95 

DIMENSION  TAB (10, 10)  96 

IF(IENT.NE.l)  GO  TO  1003  97 

CALL  RANDO ( N 1 , N 1 2 , XRN )  98 

N1=N12  99 

1 1  =  I  NT  ( 1 0 . 0*XRN )  + 1  ]  oo 

CALL  RANDO (N2,N22, XRN)  10 1 

N2=N22  1Q2 
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JJ=INT ( 10 . 0*XRN)+1 

103 

DO  1000  1=1,10 

104 

DO  1001  J=l,10 

105 

CALL  RANDO (NR , NR2 , XRN ) 

106 

NR=NR2 

107 

1001 

TAB (I , J)=XRN 

108 

1000 

CONTINUE 

109 

IC=1 

110 

1003 

DRAW=TAB (II , JJ) 

111 

CALL  RANDO (NR, NR2, XRN) 

112 

NR=NR2 

113 

TAB (II, JJ)=XRN 

114 

IF (MOD ( IC , 2) . EQ . 0)  GO  TO  1005 

115 

CALL  RANDO (N1,N 12, XRN) 

116 

N1=N12 

117 

II=INT( 10 . 0*XRN)+1 

118 

GO  TO  1006 

119 

1005 

CALL  RANDO (N2.N22, XRN) 

120 

N2=N22 

121 

J  J= I NT ( 1 0 . 0* XRN )  + 1 

jn 

1006 

IC=IC+1 

123 

RETURN 

124 

END 

125 

SUBROUTINE  RANDO ( IX , IY , YFL) 

126 

IY=IX*65539 


127 
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IF(IY)  5,6,6 

5  I Y=IY+2 147483647+1 

6  YFL=FLOAT ( IY)*0 . 4656613E-9 
RETURN 

END 

**************************  DATA 
SIMULATION  WITHOUT  ZEROES  PERIOD  1 


750367995392145 


100 


15 


.25925  .44244  1.00052  1.86506 

14.78448  11.89578  2.0 


128 

129 

130 

131 

132 

133 

134 

135 

136 

137 
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Notes  For  Program  3 

Line  # 

Comment 

19 

Initialize  the  random  number  matrix. 

20  -  24 

Compute  the  sliding  polynomial  coefficients  for  3  arcs 

25 

Range  of  v  for  one  arc. 

28  -  31 

Draw  a  random  number  and  find  which  arc  it  falls  in. 

32  -  42 

Initialize  appropriate  arc  for  interval  - ha  1 vi ng  method 

43  -  58 

Compute  a  value  H(I)  by  reverse  interpolation  from 
random  R,  using  interval  -  ha  1 ving  method. 

59  -  67 

Place  sample  items  in  order  of  magnitude. 

67  -  71 

Print  ranked  sample. 

72  -  74 

Store  ranked  sample  in  array. 

76  -  86 

Place  each  sample  rank  in  order  of  magnitude  across 
sampl es . 

87  -  92 

Print  ranked  sample  ranks. 

95  -  131 

Random  number  generator.  Random  draw  of  a  random 
number  from  a  10  x  10  matrix,  and  replacement. 

Sample  Output  For  Program 
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41.084  2  31.451  4  30.999  4  30.179  5  26.293  6  24.368  7  11.469  8  9882 

9.305  10  7.680  11  7.110  12  3.821  13  3.674  14  2.156  15  1.841 
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Sample  Output  For  Program  3  (continued) 
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PROGRAM  4 


Simulation  of  a  number  (NS)  of  synthetic  samples  of  a  number  (Nl)  of 
items  each. 


Sample  contains  zeroes. 


Nodes  were  derived  using  Program  2. 
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C  SIMULATION  WITH  ZEROES  --  USE  WITH  PROGRAM  NO .  2  1 

DIMENSION  H(IOO) ,0(3,4) ,P(6) ,S0RT(110, 20) ,AT(80)  2 

READ (5 , 8080 )  (AT(I ) , 1=1 , SO)  3 

8080  FORMAT (80A1)  4 

WRITE (6,8081)  (AT(I ) , 1=1 , 80)  5 

8081  FORMAT ('  ' ,80A1)  6 

READ (5 , 1000 )  NR ,N1 ,N2 ,NS ,NI  7 

1000  FORMAT (315,2110)  8 

WRITE (6,50)  NR,N1 ,N2 ,NS ,NI  9 

50  FORMAT ( '  ’ ,6X/ ,518)  10 

READ (5 , 1001)  (P(I ) , 1=3 , 6)  11 

1001  FORMAT (4F10 . 6 )  12 

WRITE (6,1001)  (P ( I ) ,1=3,6)  13 

READ (5 , 1001)  HB,SD,Q  14 

WRITE(6, 1001)  HB,SD,Q  15 

P ( 1 )=P ( 3 )  16 

P (2 )=0 . 0  17 

CALL  RAND (NR , N 1 , N2 , R , 1 )  18 

DO  1002  1=1,3  19 

C(I,l)=(9*(P(I+l)+P(I+2))-P(I)-P(I+3))/16  20 

C(I,2)=(ll*(P(I+2)-P(I+l))+P(I)-P(I+3))/8  21 

C ( I ,3)=(P(I)-P(I+l)-P(I+2)+P(I+3))/4  22 

1002  C(I,4)=(3*(P(I+l)-P(I+2))-P(I)+P(I+3))/2  23 

VS=(4.-2.5*EXP(-.91629*HB/Q/SD))/3.0  24 


DO  1003  M=1 ,NS 


25 


52 

DO  1004  1=1 ,NI  26 

CALL  RAND (NR, N1 ,N2 ,R, 2)  27 

IF  (R . LT . P (3 ) )  GO  TO  1005  28 

IF  (R . LT . P (4) )  GO  TO  1006  29 

IF  (R.LT.P(5))  GO  TO  1007  30 

H( I )=0 . 0  31 

GO  TO  1004  32 

1005  BL=0 . 0  33 

BR=VS  34 

IS=1  35 

GO  TO  1008  36 

1006  BL=VS  37 

BR=2 . 0*VS  38 

IS=2  39 

GO  TO  1008  40 

1007  BL=2 . 0*VS  41 

BR=3 . 0*VS  42 

IS=3  43 

1008  A=(-BL)/(BR-BL)-0.5  44 

B=l/ (BR-BL)  45 

1100  Z=A+B*(BL+BR)/2.0  46 

PC=C (C (IS , 4)*Z+C (IS , 3) )*Z+C ( IS , 2 ) )*Z+C ( IS , 1 )  47 

IF(R.LT.PC)  GO  TO  1009  48 

BL=(Z-A)/B  49 

GO  TO  1010  50 

51 


1009  BR=(Z-A)/B 
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1010  IF (ABS (BR-BL) .LT. 0,001)  GO  TO  1011 
GO  TO  1100 

1011  V=(BL+BR)/2 

IF (V . LT . 1.5)  GO  TO  1012 
H ( I ) =HB+Q*SD*ALOG ((4.0-V)/2.5)/0.91629 
GO  TO  1004 

1012  H ( I ) =HB - Q*SD*ALOG (V/1.5)/l. 52715 
1004  CONTINUE 

DO  10  J=1,NI 
X6=H(J) 

DO  20  Jl=l ,NI 

IF (H ( J1 ) . GE.X6)  GO  TO  20 

H(J)=H(J1) 

H( J1)=X6 
X6=H ( J) 

20  CONTINUE 
10  CONTINUE 

WRITE (6, 1017)  M 

1017  FORMAT ('  '/, 10X SAMPLE  NO.  ',14) 
WRITE(6, 1018)  ( K , H ( K ) , K= 1 , N I ) 

1018  FORMAT ('  ’ , 8 ( 14 , F8 . 3 ) / ) 

DO  2000  K=1,NI 

SORT(M , K)=H (K) 

2000  CONTINUE 
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57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 


75 


1003  CONTINUE 


76 
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DO  2010  K=1 ,NI  77 

DO  2005  1=1, NS  78 

CHECK=SORT(I,K)  79 

DO  2004  J=1 , NS  80 

IF(SORT(J,K) .GE. CHECK)  GO  TO  2004  81 

SORT ( I ,K)=SORT( J , K)  82 

SORT( J , K)=CHECK  83 

CHECK=SORT ( I , K)  84 

2004  CONTINUE  85 

2005  CONTINUE  86 

2010  CONTINUE  87 

WRITE ( 6 , 2223 )  (L,L=1,NI)  88 

2223  FORMAT ('  * , / / 1 IX , '  SORTED  DATA  ' / , 3X , 20 ( 12 , 4X) )  89 

DO  2001  M=1 , NS  90 

WRITE (6,2222)  (SORT(M,K) ,K=1 ,NI )  91 

2222  FORMATC  '  ,20F6.2)  92 

2001  CONTINUE  93 

70  STOP  94 

END  95 

SUBROUTINE  RAND (NR , N1 , N2 , DRAW , IENT)  96 

DIMENSION  TAB (10, 10)  97 

IF ( IENT . NE . 1 )  GO  TO  1003  98 

CALL  R ANDO ( N 1 , N 1 2 , XRN )  99 

N1=N12  100 

1 1  =  I  NT  ( 1 0 . 0  "'XRN )  + 1  101 


CALL  RANDO (N2 , N22 , XRN) 


102 
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N2=N22 

J  J= I NT ( 1 0 . 0*  XRN )  + 1 
DO  1000  1=1,10 
DO  1001  J=l,10 
CALL  RANDO ( NR , NR2 , XRN ) 

NR=NR2 

1001  TAB ( I , J ) =XRN 
1000  CONTINUE 
IC=1 

1003  DRAW=TAB (II , JJ) 

CALL  RANDO ( NR , NR2 , XRN ) 

NR=NR2 

TAB (II, JJ)=XRN 

IF (MOD ( IC , 2 ) . EQ . 0 )  GO  TO  1005 
CALL  RANDO ( N 1 , N 1 2 , XRN ) 

N1=N12 

1 1  =  I  NT  ( 1 0  .  O'*  XRN )  + 1 
GO  TO  1006 

1005  CALL  RANDO (N2.N22, XRN) 

N2=N22 

J  J= I NT ( 1 0 . O* XRN )  + 1 

1006  IC=IC+1 
RETURN 
END 
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110 
111 
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119 

120 
121 
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SUBROUTINE  RANDO ( IX , IY , YFL) 
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IY=IX*65539 
IF(IY)  5,6,6 

5  IY=IY+2 147483647+1 

6  YFL=FLOAT(IY)*0. 46566 13E-9 
RETURN 

END 

*************************  DATA  **************************** 
SIMULATION  WITH  ZEROES.  WATKINSVILLE ,  GA .  JULY 
483389174482781  100  15 

.15463  .32253  0.89981  2.81714 

2.61806  3.40995  2.0 
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Notes  For  Proqram  4 

Line  # 

Comment 

18 

Initialize  the  random  number  matrix. 

19-23 

Compute  the  Sliding  Polynomial  Coefficients. 

24 

Range  of  v  for  one  arc. 

27  -  30 

Draw  a  random  number  and  find  which  arc  it  falls  in. 

31  -  32 

Zero  if  beyond  P(5) . 

33  -  43 

Initialize  appropriate  arc  for  interval -hal ving  method 

44  -  59 

Compute  a  value  H(I)  by  reverse  interpolation  from 
random  R,  using  interval -hal vi ng  method. 

60  -  68 

Place  sample  in  order  of  magnitude. 

69  -  72 

Print  ranked  sample. 

73  -  75 

Store  ranked  sample  in  array. 

77  -  87 

Place  each  sample  rank  in  order  of  magnitude  across 
sampl es . 

88  -  92 

Print  ranked  sample  ranks. 

96  -132 

Random  number  generator.  Random  draw  of  a  random 
number  from  a  10  x  10  matrix,  and  replacement. 

Sample  Output  For  Program 
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12-903  2  11.983  3  9.868  4  9.457  5  5.062  6  2.685  7  2.553  8  2  086 

1.239  10  0.623  11  0.520  12  0.172  13  0.039  14  0.0  15  0.0 


Sample  No.  6  -  Sample  No.  95  Not  Shown 
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